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Abstract An important issue in the tomographic reconstruction of the solar poles 
is the relatively rapid evolution of the polar plumes. We demonstrate that it is 
possible to take into account this temporal evolution in the reconstruction. The 
difficulty of this problem comes from the fact that we want a 4D reconstruction (three 
spatial dimensions plus time) while we only have 3D data (2D images plus time). 
To overcome this difficulty, we introduce a model that describes polar plumes as 
stationary objects whose intensity varies homogeneously with time. This assumption 
can be physically justified if one accepts the stability of the magnetic structure. This 
model leads to a bilinear inverse problem. We describe how to extend linear inversion 
methods to these kinds of problems. Studies of simulations show the reliability of our 
method. Results for SOHO/EIT data show that we are able to estimate the temporal 
evolution of polar plumes in order to improve the reconstruction of the solar poles 
from only one point of view. We expect further improvements from STEREO/EUVI 
data when the two probes will be separated by about 60° . 



1. Introduction 

A method known as solar rotational tomography has been used to retrieve the 3D 
geometry of the solar corona (Frazin 2000; Frazin and Janzen 2002). This method 
assumes the stability of the structures during the time necessary to acquire the 
data. Since we generally have only one point of view at our disposal, about 15 
days are required to have data for half a solar rotation at the poles. Here, we focus 
our study on solar polar plumes. They are bright, radial, coronal ray structures 
located at the solar poles in regions of open magnetic field. The study of plumes is 
of great interest since it may be the key to the understanding of the acceleration 
of the fast component of the solar wind (Teriaca et al, 2003). However the three- 
dimensional shape of these structures is poorly known and different assumptions have 
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been made, e.g. Gabriel et al, 2005; Llebaria, Saez, and Lamy, 2002. The plumes 
are known to evolve with a characteristic time of approximately 24 hours on spatial 
scales typical of Extreme ultra-violet Imaging Telescope (SOHO/EIT) data (2400 
km) (DeForest, Lamy, and Llebaria, 2001). Consequently the stability assumption 
made in rotational tomography fails. Fortunately, the Solar TErestrial RElations 
Observatory (STEREO) mission consists of two identical spacecraft STEREOa and 
STEREOb which take pictures of the Sun from two different points of view. With the 
SOHO mission still operating, this results in three, simultaneous points of view. Three 
viewpoints help to improve the reconstruction of the plumes, but they are still not 
enough to use standard tomographic algorithms. The problem is underdetermined 
and consequently one has to add a priori information in order to overcome the lack 
of information. This leads to challenging and innovative signal analysis problems. 
There are different ways to deal with underdetermination depending on the kind of 
object to be reconstructed. Interestingly the field of medical imaging faces the same 
kind of issues. In cardiac reconstruction, authors make use of the motion periodicity 
in association with a high redundancy of the data (Grass et al, 2003; Kachelriess, 
Ulzheimer, and Kalender, 2000). If one can model the motion as an afRne transfor- 
mation, and if one assumes that we know this transformation, one can obtain an 
analytic solution (Ritchie et al, 1996; Roux et al, 2004). 

In solar tomography, the proposed innovative approaches involve the use of ad- 
ditional data such as magnetic-field measurements in the photosphere (Wiegelmann 
and Inhester, 2003) or data fusion (Frazin and Kamalabadi, 2005). Attempts have 
been made by Frazin et al. (2005) to treat temporal evolution using Kalman filtering. 

Since polar plumes have apparently a local, rapid, and aperiodic temporal evo- 
lution, we developed as in the previously referenced work, a model based on the 
specifics of the object we intend to reconstruct (preliminary results can be found in 
Barbey et al, (2007). Plumes have an intensity which evolves rapidly with time, 
but their position can be considered as constant. This hypothesis is confirmed by 
previous studies of the plumes such as DeForest, Lamy, and Llebaria (2001). The 
model is made up of an invariant morphological part (x) multiplied by a gain term 
{Ot) that varies with time. Only one gain term is associated with each plume in order 
to constrain the model. So we assume that the position of each plume in the scene 
is known. This model is justified if we consider polar plumes to be slowly evolving 
magnetic structures in which plasma flows. 

Thanks to this model we can perform time-evolving three-dimensional tomogra- 
phy of the solar corona using only extreme ultra-violet images. Furthermore, there is 
no complex, underlying physical model. The only assumptions are the smoothness of 
the solution, the area-dependant evolution model, and the knowledge of the plume 
position. These assumptions allow us to consider a temporal variation of a few days, 
while assuming only temporal smoothness would limit variations to the order of one 
solar rotation (about 27 days). To our knowledge, the estimation of the temporal 
evolution has never been undertaken in tomographic reconstruction of the solar 
corona. 

We first explain our reconstruction method in a Bayesian framework (Section 
2). We then test the validity of our algorithm with simulated data (Section 3). An 
example of a reconstruction on real SOHO/EIT data is shown in Section 4. Results 
are discussed in Section 5. We conclude in Section 6. 
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Figure 1. Scheme of the data acquisition geometry. (O; x,y, z) defines the Carrington heUo- 
centric frame of reference. S is the spacecraft considered. (j> is the latitude, and 6 the longitude 
of this spacecraft. V is the virtual detector. 

2. Method 

Tomographic reconstruction can be seen as an inverse problem, the direct problem 
being the acquisition of data images knowing the emission volume density of the 
object (Section 2.1). If the object is evolving during the data acquisition, the inverse 
problem is highly underdetermined. So our first step is to redefine the direct problem 
thanks to a reparametrization, in order to be able to define more constraints (Section 
2.2). Then, we place ourselves in the Bayesian inference framework in which data 
and unknowns are considered to be random variables. The solution of the inverse 
problem is chosen to be the maximum a posteriori (Section 2.3). This leads to a 
criterion that we minimize with an alternate optimization algorithm (Section 2.4). 

2.1. Direct Problem 

The geometrical acquisition is mathematically equivalent to a conical beam data 
acquisition with a virtual spherical detector (see Figure 1). In other words, the step 
between two pixels vertically and horizontally is constant in angle. The angle of the 
full field of view is around 45 minutes. In order to obtain an accurate reconstruc- 
tion, we take into account the exact geometry, which means the exact position and 
orientation of the spacecraft relatively to Sun center. We approximate integration of 
the emission in a fiux tube related to a pixel by an integration along the line of sight 
going through the middle of that pixel. We choose to discretize the object in the 
usual cubic voxels, a; is a vector of size A'' containing the values of all voxels. In the 
same way, we define the vector of data yi, of size M at time t. Since the integration 
operator is linear, the projection can be described by a matrix Pt. We choose nt to 
be an additive noise: 



Pt is the projection matrix at time t of size M x N which is defined by the position 
and the orientation of the spacecraft at this time. Its transpose is the backprojec- 
tion matrix. Note that a uniform sampling in time is not required. In order to be 
able to handle large problems with numerous well-resolved data images and a large 
reconstruction cube, we chose not to store the whole projection matrix. Instead, 



yt = PtXt + nt,'ite [1,...,T] 



(1) 
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we perform the projection operation (Pa;) or its transpose each time it is needed 
at each iteration. Thus, we need a very efficient algorithm. We developed a code 
written in C which performs the projection operation. It makes use of the geometrical 
parameters given in the data headers in order to take into account the exact geometry 
(conicity, position, and orientation of the spacecraft). To keep this operation fast, 
we implemented the Siddon algorithm (Siddon, 1985). It allows a fast projection or 
backprojection in the case of cubic voxels (Cartesian grid). Since we focus on a small 
region at the poles, we consider that we do not need to use a spherical grid which 
would require a more time-consuming projection algorithm. 

We take into account the fact that the field of view is conical. Despite the fact 
that the acquisition is very close to the parallel acquisition geometry, it is sufficient 
to introduce an error of several voxels of size 0.01 solar radius from one side to the 
other of a three solar radii reconstructed cube. 

2.2. Modeling of the Temporal Evolution 

With this model, the inverse problem is underdetermined since we have at most 
three images at one time and we want to reconstruct the object with its temporal 
evolution. In order to do so, we first redefine our unknowns to separate temporal 
evolution from spatial structure. We introduce a new set of variables gt of size A'' 
describing the temporal evolution and require that x does not depend on time: 



with o being the term- by-term multiplication of vectors. This operator is clearly 
bilinear. However, this model would increase the number of variables excessively. 

So, we need to introduce some other kind of a priori into our model. We make 
the hypothesis that all of the voxels of one polar plume have the same temporal 
evolution: 



The matrix L of size N x P (P being the number of areas) localizes areas where 
the temporal evolution is identical. Each column of L is the support function of one 
of the plumes. We would like to stress that in our hypothesis, those areas do not 
move relative to the object. In other words, L does not depend on time. Localizing 
these areas defines L and only leaves PT variables to estimate. We redefined our 
problem in a way that limits the number of parameters to estimate but still allows 
many solutions. Furthermore, the problem is linear in x knowing and linear in 
knowing x. It will simplify the inversion of the problem as we shall see later. Note, 
however that the uniqueness of a solution (x, 6) is not guaranteed with bilinearity 
despite its being guaranteed in the linear case. This example shows that A can be 
chosen arbitrarily without changing the closeness to the data: xog = (Ax) o (^A~^g) , 
where A is a real constant. Introducing an a priori of closeness to 1 for would allow 
us to deal with this indeterminacy in principle. But note that this indeterminacy is 
not critical since the physical quantity of interest is only the product xog. Feron, 
Duchene, and Mohammad-Djafari (2005) present a method which solves a bilinear 
inversion problem in the context of microwave tomography. 

We do not deal with the estimation of the areas undergoing evolution, but we 
assume in this paper that the localization is known. This localization can be achieved 



yt = Pt{x o gt) + nt 



(2) 



gt = LOt 



(3) 
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using other sources of information, e.g. stereoscopic observations. We expect to be 
able to locate the areas using some other source of information. 

We can regroup the equations of the direct problem. We have two ways to do so, 
each emphasizing the linearity throughout one set of variables. 



"1 


■( 







y = UxO + n 
PiXL \ ( Oi 



PtXl ) \eT 



+ : 



(4) 



with X = diag{x), the diagonal matrix defined by x. x is of size N, y and n are of 
size MT,0 is of size PT and Ux is of size MTx PT. 
Similarly, 



y = Vex + n 
Pidie.g{L0i) 



with Vg 




(5) 



PTdiag(L6>T) J \ Id 

with Id the identity matrix of size M x M. Ve is of size MT x TV. 

2.3. Inverse Problem 

In Bayes' formalism, solving an inverse problem consists in knowing the a posteriori 
(the conditional probability density function of the parameters, the data being given). 
To do so we need to know the likelihood (the conditional probability density function 
of the data knowing the parameters) and the a priori (the probability density func- 
tion of the parameters). An appropriate model is a Gaussian, independent, identically 
distributed (with the same variance) noise n. The likelihood function is deduced from 
the noise statistic: 

f{y\x, e, an,M) = exp ) (6) 

Ai = [P, L] describing our model (the projection algorithm and parameters and the 
choice of the plume position). We assume that the solution is smooth spatially and 
temporally, so we write the a priori as follows: 

f{x\ax) = K2e^p (-^^3^) = ^3exp ("^f^) (7) 

Dr and Z?t axe discrete differential operators in space and time. Bayes' theorem 
gives us the a posteriori law if we assume that the model Ai is known as well as the 
hyperparameters H = [an, o^x, ere]' 

f{x,e\y,n,M)- fiy\n,M) 
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We need to choose an estimator. It allows us to define a unique solution instead of 
having a whole probability density function. We then choose to define our solution 
as the maximum a posteriori, which is given by: 

(x"*'',^"*'") = a.rg max f{y\x,e, an, M)f{x\ax)f{0\cTe) (9) 

x,e 

since f{y\A4) is a constant. Equation (9) can be rewritten as a minimization problem: 

^^MAP^^MAP^ = argmin J(x,6>) (10) 

x,0 



with: 



J{x,e) = -2anlogfix,e\y,M,n) = Wv-UxOf + X\\Drxf +n\\Dtef (11) 

A = ^ and n= are user-defined hyperparameters. 

The equivalence of Equations (9) and (10) has been proved by Demoment (1989). 

Note that the solution does not have to be very smooth. It mostly depends on 
the level of noise since noise increases the underdetermination of the problem as it 
has been shown by the definition of A and ji. 



2.4. Criterion Minimization 



The two sets of variables x and 6 are very different in nature. However, thanks 
to the problem's bilinearity, one can easily estimate one set while the other is fixed. 
Consequently we perform an iterative minimization of the criterion, and we alternate 
minimization of x and 6. At each step n we perform: 

6>"+^ = argmin J(x", 6)) and = argmin J (x, 0"+^) (12) 

X 

The two subproblems are formally identical. However, 6 is much smaller than x. 

This is of the utmost practical importance since one can directly find the solution on 
9 by using the pseudo-inverse method, x is too big for this method, and we have to 
use an iterative scheme such as the conjugate-gradient to approximate the minimum. 
These standard methods are detailed in Appendices A and B. 



2.5. Descent Direction Definition and Stop Threshold 



We choose to use an approximation of the conjugate-gradient method that is known 
to converge much more rapidly than the simple gradient method (Nocedal and 
Wright, 2000; Polak and Ribiere, 1969). 

MP _ (v^JL=.P,V;.JL=^p-i) (13) 

^ - l|V.J|_p-i|P 

Since the minimum is only approximately found, we need to define a threshold 
which we consider to correspond to an appropriate closeness to the data in order to 
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stop the iterations. Since the solution is the point at which the gradient is zero, we 
choose this threshold for updating x: 

mean^g[a,p ,^p-i^p-2]||VxJ||^ < Sx (14) 

For the global minimization, the gradient is not computed, so we choose: 

meanjjj ,j_2] || (xji, en)-{Xn-l,en-l)r <Sg (15) 

Note that this way to stop the iteration allows one to define how close one wants to 
be to the solution: if the difference between two steps is below this threshold, it is 
considered negligible. The algorithm can be summarized as shown in Figure 2. 

initialize : a; = and = 1 
while Equation (15) is satisfied 

X minimization: 

while Equation (14) is satisfied 

* compute gradient at Xn with Equation (20) 

* compute descent direction with Equation (13) 

* compute optimum step with Equation (22) 

* update X with Equation (23) 
endwhile 

minimization: 

* compute the matrix U^nUx" and the vector U^ny 

* inverse the matrix Ux'* + ^iDJ Dr 

* compute Equation (19) 

endwhile 

Figure 2. Tomographic Reconstruction with Temporal Evolution Algorithm 



3. Method Validation 

In order to validate the principle of our method and test its Umits, we simulate an 
object containing some plumes with temporal evolution and try to extract it from 
the data. 

3.1. Simulation Generation Process 

We generate an emission cube with randomly-placed, ellipsoidal plumes with a 
Gaussian shape along each axis: 



Ep = A exp 



1 


r.u^) 


' 1 




2 


a 


2 


V b \ 



(16) 



The plumes evolve randomly but smoothly by interpolating over a few randomly 
generated points. Once the object is generated, we compute a typical set of 60 images 
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equally spaced along 180° using our projector algorithm. A Gaussian random noise 
is added to the projections with a signal to noise ratio (SNR) of five. The simulation 
parameters are summarized in Table 1. 



Table 1. Simulation Definition: Plumes Pcirameters 



Plume 
Number 


Scmimajor 
Axis a 


Semiminor 
Axis b 


4> 


Xo 


yo 


Intensity 
(A) 


1 


4.8 


4.2 


1.2 


29 


29 


329 


2 


5.6 


3.3 


1.1 


23 


33 


430 


3 


5.2 


4.8 


0.1 


40 


42 


723 



3.2. Results Analysis 

We now compare our results (Figure 3) with a filtered back-projection (FBP) algo- 
rithm. This method is explained by Natterer (1986) and Kak and Slaney (1987). 

By comparing the simulation and the reconstruction in Figure 3, we can see the 
quality of the temporal evolution estimation. The shape of the intensity curves is well 
reproduced except for the first plume in the first ten time steps where the intensitj^ is 
slightly underestimated. This corresponds to a period when plume 1 is hidden behind 
plume 2. Thus, our algorithm attributes part of the plume 1 intensity to plume 2. 
Let us note that this kind of ambiguity will not arise in the case of observations from 
multiple points of view such as STEREO/EUVI observations. The indeterminacy of 
the problem is due to its bilinearity discussed in Section 2.2. This allows the algorithm 
to attribute larger values to the 6 parameters and to compensate by decreasing the 
corresponding x. This is not a drawback of the method since it allows discontinuities 
between plumes and interplumes. The only physical value of interest is the product 

« Off- 
Figure 4 shows the relative intensity of the plumes at different times. One can 
compare with the reconstruction. One way to quantify the quality of the reconstruc- 
tion is to compute the distance (quadratic norm of the difference) between the real 
object and the reconstructed one. Since the FBP reconstruction does not actually 
correspond to a reconstruction at one time, we evaluate the minimum of the distances 
at each time. We find it to be 3000. This is to be compared with a value of 700 with 
our algorithm, which is much better. 



Table 2. Simulation Definition: Geometric Parameters 



cube size 


cube number 


pixel 


projection 


(solar radii) 


of voxels 


size (radians) 


number of pixels 


1 X 1 X 0.05 


64 X 64 X 4 


5 X 10-5 X 5 X 10-5 


128 X 8 
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X axis in solar radii 



X axis in solar radii 



X axis in solar radii 



(a) X with a FBP algorithm 



(b) Simulated x 




X axis in solar radii 




(c) X with our algorithm 




(d) Chosen temporal areas 



(e) Simulated 



(f) 6 with our algorithm 



Figure 3. Comparison of a standard FBP method (a), the real simulated object (b), and the 
object reconstructed with our method (c). The object is reconstructed using 60 projections 
regularly spaced over 180° . The areas of homogeneous temporal evolution (e) are the same in 
the simulation and the reconstruction. We associated one time per projection to define in 
the simulation (e) and our reconstruction (f). The time scale is in days assuming a rotation 
speed of half a rotation in 14 days, x is the spatial distribution of the emission density volume. 
is a gain representing the emission variation over time. Except for the FBP reconstruction, 
only the product x o has physical dimensions. The spatial scales are given in solar radii 
and centered on the solar axis of rotation, (a), (b) and (c) are slices of 3D cubes at the same 
z = 0.1 Rq. Emission densities (arbitrary units) are scaled in the color bars in the right-end 
side of (a), (b), (c). 




X axis in solar radii 



X axis in solar radi 



X axis in solar radii 



(a) Simulation at lOAT 

-0.5| 



(b) Simulation at 30AT 



(c) Simulation at 50AT 




X axis in solar radi 



(d) Result at lOAT 



(e) Result at 30AT 



(f) Result at 50AT 



Figure 4. Comparison oi xog simulated and reconstructed at different times. AT is the time 
between two data images (5.6 hours). Distances are in solar radii. Values represent the volume 
emission density. All of this images are slices of 3D cubes at the same z = O.IRq. 
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Table 3. Simulation Definition: Other Parameters 



SNR 


A 




5'x 


Sg 


5 


2 X IQ-^ 


100 


2 X IQ-^ 


1 X 10^2 




(a) New choice of areas (b) x with our algorithm (c) with our algorithm 

Figure 5. Reconstruction with smaller areas. To be compared with Figure 3. The new areas 
(a) do not correspond anymore to the ones used to generate the data, (b) is the emission map 
and (c) the temporal evolution estimated with our algorithm, (b) and (c) are slices of 3D cubes 
at the same z = 0.1 Rq. Emission densities (arbitrary units) are scaled in the color bars in the 
right-end side of (b). 



3.3. Choice of Evolution Areas 

One can think that the choice of the evolution areas is critical to the good perfor- 
mance of our method. We show in this section that it is not necessarily the case by 
performing a reconstruction based on simulations with incorrect evolution areas. All 
parameters and data are exactly the same as in the previous reconstruction. The 
only difference is in the choice of the areas, i.e. the L matrix. These are now defined 
as shown in Figure 5(a). 

Although approximately 50 % of the voxels are not associated with their correct 
area, we can observe that the algorithm still performs well. The emission map of 
Figure 5(b) is still better than the emission reconstructed by a FBP method. Plus, 
the estimation of the temporal evolution in Figure 5(c) corresponds to the true 
evolution 3(e) even if less precisely than in Figure 3(f). 



4. Reconstruction of SOHO/EIT Data 

4.1. Data Preprocessing 

We now perform reconstruction using SOHO /EIT data. We have to be careful when 
applying our algorithm to real data. Some problems may arise due to phenomena 
not taken into account in our model; e.g. cosmic rays, or missing data. 

Some of these problems can be handled with simple preprocessing. We consider 
pixels hit by cosmic rays as missing data. They are detected with a median filter. 
These pixels and missing blocks are labeled as missing data and the projector and 
the backprojector do not take them into account {i.e. the corresponding rows in the 
matrices are removed). 
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Table 4. EIT Data Reconstruction: Geometric Parameters 



cube size 


cube number 


pixel 


projection 


(solcir radii) 


of voxels 


size (radians) 


number of pixels 


3 X 3 X 0.15 


256 X 256 X 8 


2.55 X 10-5 X 2.55 X 10-5 


512 X 38 



Table 5. EIT Data Reconstruction 
: Other Parameters 



A 


M 




So 


2 X 10-2 


1 X 10'' 


0.1 


0.05 



4.2. Results Analysis 

In Figures 6 and 7, we present results from 17.1 nm EIT data between 1 and 14 
November 1996. This period corresponds to the minimum of solar activity when 
one can expect to have less temporal evolution. 17.1 nm is the wavelength where 
the contrast of the plumes is the strongest. Some images are removed resulting in a 
sequence of 57 irregularly-spaced projections for a total coverage of 191°. We assume 
that we know the position of four evolving plumes as shown on Figure &(h). For each 
reconstructed image, we present subareas of the reconstructed cube of size 64x64 
centered on the axis of rotation. We assume the rotation speed to be the rigid body 
Carrington rotation. All of the parameters given in Table 4 and 5 are shared by the 
different algorithms provided they are required by the method. The computation of 
this reconstruction on a Intel(R) Pentium(R) 4 CPU 3.00 GHz was 13.5 hours long. 

Presence of negative values is the indication of a poor behavior of the tomographic 
algorithm since it does not correspond to actual physical values. We can see in Figure 
6 that our reconstruction has many fewer negative values in the x map than the FBP 
reconstruction. In the FBP reconstruction cube, 50% of the voxels have negative 
values; in the gradient-like reconstruction without temporal evolution 36% of the 
voxels are negative while in our reconstruction only 25 % are negative. This still 
seems like a lot but most of these voxels are in the outer part of the reconstructed 
cube. The average value of the negative voxels is much smaller also. It is -120 for 
the FBP, -52 for the gradient-like method without temporal evolution, and only -19 
for our reconstruction with temporal evolution. However, we notice that the gain 
coefficients present a few slightly negative values. 

In the reconstructions without temporal evolution, plumes three (upper right) 
and four (lower right) correspond to a unique elongated structure which we choose 
to divide. Note how our algorithm updated the x map reducing the emission values 
between these two plumes. It shows that what was seen as a unique structure was an 
artifact resulting from temporal evolution and it tends to validate the usefulness of 
our model. We note the disappearance of a plume located around (-0.2, -0.15) solar 
radii on the FBP reconstruction. It shows the utility of gradient-like methods to get 
rid of artifacts due to the non-uniform distribution of images. Another plume at (0.2, 
0.2) solar radii has more intensity in the reconstruction without temporal evolution 
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Figure 6. A comparison of FBP (a), a gradient-like algorithm without temporal evolution 
(c), and our algorithm (d) with real EIT data, x is the spatial distribution of the volume 
emission density integrated over EIT 17.1 nm passband. The chosen areas are shown in (b). 
is a gain representing the emission variation during time (e). The time scale is in days. In 
the case of our algorithm, only the product x o0 has physical meaning. The spatial scales are 
given in solar radii and centered on the solar axis of rotation, (a), (b), (c), and (d) are slices of 
3D cubes at the same z = 1.3-Rq. Emission densities (arbitrary units) are scaled in the color 
bars in the right-end side of (a), (c), (d). 
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(a) Reconstruction at i25 (b) Reconstruction at {35 (c) Reconstruction at ti^ 
(6 days) (8.5 days) (9.8 days) 

Figure 7. Reconstruction o{ x o g at different times. Distances are in solar radii. Values 
represent the volume emission density integrated over the EIT 17.1 nm passband. All of these 
images are slices of 3D cubes at the same ^ = 1.3 Rq. 

than with our algorithm. It illustrates how temporal evolution can influence the 
spatial reconstruction. 

5. Discussion 

The major feature of our approach is the quality of our reconstruction, which is 
much improved with respect to FBP reconstruction, as demonstrated by the smaller 
number of negative values and the increased closeness to the data. Let us now discuss 
the various assumptions that have been made through the different steps of the 
method. 

The strongest assumption we made, in order to estimate the temporal evolution of 
polar plumes, is the knowledge of the plume position. Here, we choose to define the 
plumes as being the brightest points in a reconstruction without temporal evolution. 
The choice is not based on any kind of automatic threshold. The areas are entirely 
hand-chosen by looking at a reconstruction. It is possible that these areas do not 
correspond to the actual physical plumes, they could correspond to areas presenting 
increased emission during half a rotation. Note that this is biased in favor of plumes 
closer to the axis of rotation since, along one slice of the reconstructed cartesian 
cube, their altitude is lower and thus, their intensity is higher. In order to have 
constant altitude maps one would have to carry out the computation on a spherical 
grid or to interpolate afterwards onto such a grid. For this reconstruction example 
we are aware that we did not locate all of the plumes but only tried to find a few. It 
would be interesting to try to locate the plumes using other data or with a method 
estimating their positions and shapes. 

The method involves hyperparameters which we choose to set manually. There 
are methods to estimate hyperparameters automatically such as the L-curve method, 
the cross-validation method (Golub, Heath, and Wahba, 1979) or the fuU-bayesian 
method (Higdon et al, 1997; Champagnat, Goussard, and Idler, 1996). We per- 
formed reconstructions using different hyperparameter values. We then looked at 
the reconstruction to see if the smoothness seemed exaggerated or if the noise were 
amplified in the results. This allowed us to reduce the computational cost and does 
not really put the validity of the method into question. 

One possible issue with this algorithm is the non-convexity of our criterion. This 
can lead to the convergence to a local minimum that does not correspond to the 
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desired solution defined as the global minimum of the criterion. One way to test this 
would be to change the initialization many times. 

We chose the speed of rotation of the poles to be the Carrington rotation speed. 
But the speed of the polar structures has not been measured precisely to our knowl- 
edge and could affect drastically the reconstruction. This is an issue shared by all 
tomographic reconstructions of the Sun. 

In the current approach, we need to choose on our own the position of the time- 
evolving areas which are assumed to be plumes. This is done by assuming that more 
intense areas of a reconstruction without temporal evolution correspond to plume 
positions. A more rigorous way would be to try to use other sources of information to 
try to localize the plumes. Another, self-consistent way, would be to develop a method 
that jointly estimates the position of the plumes in addition to the emission (x) and 
the time evolution (0). We could try to use the results of Yu and Fessler (2002) 
who propose an original approach in order to reconstruct a piece-wise homogeneous 
object while preserving edges. The minimization is alternated between an intensity 
map and boundary curves. The estimation of the boundary curves is made using 
level sets technics ((Yu and Fessler, 2002) and references therein). It would also be 
possible to use a Gaussian mixture model (Snoussi and Mohammad- Djafari, 2007). 

6. Conclusion 

We have described a method that takes into account the temporal evolution of polar 
plumes for tomographic reconstruction near the solar poles. A simple reconstruction 
based on simulations demonstrates the feasibility of the method and its efficiency in 
estimating the temporal evolution assuming that parameters such as plume position 
or rotation speed are known. Finally we show that it is possible to estimate the 
temporal evolution of the polar plumes with real data. 

In this study we limited ourselves to reconstruction of images at 17.1 nm but 
one can perform reconstructions at 19.5 nm and 28.4 nm as well. It would allow us 
to estimate the temperatures of the electrons as in Frazin, Kamalabadi, and Weber 
(2005) or Barbey et al. (2006). 
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Appendix 

A. Pseudo-Inverse Minimization 

We want to minimize: 



The second term does not depend on 9. Due to the strict convexity of the crite- 
rion, the solution is a zero of the gradient. Since the criterion is quadratic, one can 
explicitly determine the solution: 



J=\\y- Ua^nOf + A||r>^a;"||^ + fiWDtOf 



(17) 




(18) 
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from which we conclude: 



rjn+l 



(19) 



B. Gradient-like Method 

In this method we try to find an approximation of the minimum by decreasing the 
criterion iteratively. The problem is divided in two subproblems: searching for the 
direction and searching for the step of the descent. In gradient-like methods, the 
convergence is generally guaranteed ultimately to a local minimum. But since the 
criterion is convex, the minimum is global. To iterate, we start at an arbitrary point 
(a;*^) and go along a direction related to the gradient. The gradient at the p*'* step 
is: 

JL=^. = 2VeI+i [Vg^+.x^ -y)+ 2\DjDrxP (20) 

Once the direction is chosen, searching for the optimum step is a linear minimization 
problem of one variable: 

cioPT = argmin J(cc^ + adF~^^) (21) 

a 

which is solved by: 

2||ye„+idP+i||2 + A||D^dP+i||2 ^^^> 

We can write the iteration: 

xf+i =xP + al+TdF^^ (23) 
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